Home > sgwt_toolbox > demo > sgwt_demo1.m

sgwt_demo1

PURPOSE ^

sgwt_demo1 : SGWT for swiss roll data set

SYNOPSIS ^

function sgwt_demo1

DESCRIPTION ^

 sgwt_demo1 : SGWT for swiss roll data set

 This demo builds the SGWT for the swiss roll synthetic data set. It
 computes a set of scales adapted to the computed upper bound on the
 spectrum of the graph Laplacian, and displays the scaling function and
 the scaled wavlet kernels, as well as the corresponding frame bounds. It
 then computes the wavelets centered on a single vertex, and displays
 them. This essentally reproduces figure 3 from 
 Hammond,Vangergheynst, Gribonval 2010.

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

0001 % sgwt_demo1 : SGWT for swiss roll data set
0002 %
0003 % This demo builds the SGWT for the swiss roll synthetic data set. It
0004 % computes a set of scales adapted to the computed upper bound on the
0005 % spectrum of the graph Laplacian, and displays the scaling function and
0006 % the scaled wavlet kernels, as well as the corresponding frame bounds. It
0007 % then computes the wavelets centered on a single vertex, and displays
0008 % them. This essentally reproduces figure 3 from
0009 % Hammond,Vangergheynst, Gribonval 2010.
0010 
0011 % This file is part of the SGWT toolbox (Spectral Graph Wavelet Transform toolbox)
0012 % Copyright (C) 2010, David K. Hammond.
0013 %
0014 % The SGWT toolbox is free software: you can redistribute it and/or modify
0015 % it under the terms of the GNU General Public License as published by
0016 % the Free Software Foundation, either version 3 of the License, or
0017 % (at your option) any later version.
0018 %
0019 % The SGWT toolbox is distributed in the hope that it will be useful,
0020 % but WITHOUT ANY WARRANTY; without even the implied warranty of
0021 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
0022 % GNU General Public License for more details.
0023 %
0024 % You should have received a copy of the GNU General Public License
0025 % along with the SGWT toolbox.  If not, see <http://www.gnu.org/licenses/>.
0026 
0027 % Create swiss roll point cloud
0028 function sgwt_demo1 
0029 close all
0030 fprintf('Welcome to SGWT demo #1\n');
0031 
0032 npoints=500;
0033 fprintf('Creating Swiss Roll point cloud with %g points\n',npoints);
0034 dataparams=struct('n',npoints,'dataset',-1','noise',0,'state',0);
0035 r=create_synthetic_dataset(dataparams);
0036 x=rescale_center(r.x);
0037 
0038 
0039 fprintf('Computing edge weights and graph Laplacian\n');
0040 % Compute Weighted graph adjacency matrix, and graph Laplacian
0041 d=distanz(x);
0042 s=.1;
0043 A=exp(-d.^2/(2*s^2)); 
0044 L=full(sgwt_laplacian(A));
0045 
0046 fprintf('Measuring largest eigenvalue, lmax = ');
0047 
0048 %% Design filters for transform
0049 Nscales=4;
0050 lmax=sgwt_rough_lmax(L);
0051 fprintf('%g\n',lmax);
0052 
0053 fprintf('Designing transform in spectral domain\n'); 
0054 [g,gp,t]=sgwt_filter_design(lmax,Nscales);
0055 arange=[0 lmax];
0056 %% Display filter design in spectral domain
0057 figure
0058 sgwt_view_design(g,t,arange);
0059 ylim([0 3])
0060 set(gcf,'position',[0 780,600,300])
0061 %% Chebyshev polynomial approximation
0062 m=50; % order of polynomial approximation
0063 fprintf('Computing Chebyshev polynomials of order %g for fast transform \n',m);
0064 for k=1:numel(g)
0065   c{k}=sgwt_cheby_coeff(g{k},m,m+1,arange);
0066 end
0067 
0068 %% compute transform of delta at one vertex
0069 jcenter=32; % vertex to center wavelets to be shown
0070 fprintf('Computing forward transform of delta at vertex %g\n',jcenter);
0071 N=size(L,1);
0072 d=sgwt_delta(N,jcenter);
0073 % forward transform, using chebyshev approximation
0074 wpall=sgwt_cheby_op(d,L,c,arange);
0075 
0076 fprintf('Displaying wavelets\n');
0077 msize=100;
0078 cp=[-1.4,-16.9,3.4]; % camera position
0079 %% Visualize result
0080 
0081 % show original point
0082 ws=300;
0083 figure;
0084 xp=0; yp=ws+100;
0085 set(gcf,'position',[xp,yp,ws-10,ws+10]);
0086 scatter3(x(1,:),x(2,:),x(3,:),msize,d,'.');
0087 set(gcf,'Colormap',[.5 .5 .5;1 0 0]);
0088 clean_axes(cp);
0089 title(sprintf('Vertex %g',jcenter));
0090 
0091 % show wavelets
0092 for n=1:Nscales+1
0093     wp=wpall{n};
0094     figure
0095     xp=mod(n,3)*(ws+10);
0096     yp=(1-floor((n)/3))*(ws+100);
0097     set(gcf,'position',[xp,yp,ws-10,ws+10]);
0098     scatter3(x(1,:),x(2,:),x(3,:),msize,wp,'.');
0099     colormap jet
0100     caxis([-1 1]*max(abs(wp)));
0101     clean_axes(cp);
0102 
0103     hcb=colorbar('location','north');
0104     cxt=get(hcb,'Xtick');
0105     cxt=[cxt(1),0,cxt(end)];
0106     set(hcb,'Xtick',cxt);
0107     cpos=get(hcb,'Position');
0108     cpos(4)=.02; % make colorbar thinner
0109     set(hcb,'Position',cpos);
0110     set(hcb,'Position',[.25 .91 .6 .02]);
0111     
0112     if n==1
0113       title('Scaling function');
0114     else      
0115       title(sprintf('Wavelet at scale j=%g, t_j = %0.2f',n-1,t(end+1-(n-1))));
0116 
0117     end
0118 end
0119 
0120 
0121 function clean_axes(cp)
0122 xlim([-1 1]);ylim([-1 1]);zlim([-1 1]);
0123 set(gca,'Xtick',[-1 0 1]);
0124 set(gca,'Ytick',[-1 0 1]);
0125 set(gca,'Ztick',[-1 0 1]);
0126 axis square
0127 set(gca,'CameraPosition',cp);
0128 
0129 % rescale_center
0130 % center input data at origin, then rescale so that all coordinates
0131 % are between -1 and 1
0132 %
0133 % x should be dxN
0134 function r=rescale_center(x)
0135 N=size(x,2);
0136 d=size(x,1);
0137 x=x-repmat(mean(x,2),[1,N]);
0138 c=max(abs(x(:)));
0139 r=x/c;

Generated on Wed 13-Oct-2010 13:36:39 by m2html © 2003